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Abstract We report the recent progress on the determination of three-nucleon forces 
(3NF) in lattice QCD. We utilize the Nambu-Bethe-Salpeter (NBS) wave function to 
define the potential in quantum field theory, and extract two-nucleon forces (2NF) 
and 3NF on equal footing. The enormous computational cost for calculating multi- 
baryon correlators on the lattice is drastically reduced by developing a novel contraction 
algorithm (the unified contraction algorithm). Quantum numbers of the three-nucleon 
(3N) system are chosen to be (/, J P ) = (1/2, l/2 + ) (the triton channel), and we extract 
3NF in which three nucleons are aligned linearly with an equal spacing. Lattice QCD 
simulations are performed using Nf = 2 dynamical clover fermion configurations at 
the lattice spacing of a — 0.156 fm on a 16 3 x 32 lattice with a large quark mass 
corresponding to m-n = 1.13 GeV. Repulsive 3NF is found at short distance. 
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1 Introduction 

Nuclear (and baryonic) forces are fundamental quantities in nuclear physics, and it is 
of crucial importance to determine them from the underlying theory, Quantum Chro- 
modynamics (QCD). In particular, the determination of few-baryon forces has a huge 
impact on not only nuclear physics but also astrophysics. 

There are various phenomena where three-nucleon forces (3NF) are considered to 
play an important role, e.g., the binding energies of light nuclei pQ, the properties of 
neutron-rich nuclei and the supernova nucleosynthesis [2] and the nuclear equation 
of state (EoS) as well as the saturation point of nuclear matter [3]- Deuteron-proton 
elastic scattering experiments have also shown a clear indication of 3NF [4][5]. Recent 
observational data on the maximum mass of neutron stars [5] triggered renewed interest 
in the nuclear EoS at high density, and universal short-range repulsion for three baryons 
(nucleons and hyperons) may be needed in neutron stars with hyperon core 7,8. 
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Despite of its phenomenological importance, microscopic understanding of 3NF is 
still limited. Pioneered by Fujita and Miyazawa [9j, the long range part of 3NF has been 
commonly modeled by the two-pion exchange (27rE), particularly with the /^-resonance. 
This 2-7TE-3NF is known to have an attractive nature at long distance, and is considered 
to play a dominant role to resolve the problem that light nuclei are underbound with 
2NF only. Unfortunately, the 27rE-3NF component alone cannot represent the whole 
properties of 3NF: For instance, it would provide too strong attraction, and thus incurs 
the overbound problem for light nuclei. In order to handle this issue, an additional 
repulsive component of 3NF at short distance is often introduced. In Tucson-Melbourne 
model [10], a phenomenological cut-off parameter, A ~ 0.7 GeV, is introduced in the 
(monopole) form factor of nNN vertex [llj . Compared to A ~ 1.3 GeV obtained from 
NN scattering data |12| . one observes that short range repulsion is implicitly added 
in 3NF. In Urbana/Illinois models, on the other hand, a repulsive 3NF component at 
short distance is explicitly introduced in a purely phenomenological way [13] . Recently, 
an approach based on the chiral effective field theory (^EFT) is also developing [141115] , 
and found to be useful to classify and parametrize the two-, three- and more-nucleon 
forces. The xEFT-3NF contains not only pion exchange terms but also short range 
terms. The LECs (c£>, eg) corresponding to the latter are known to be sensitive to the 
cut-off parameter in xEFT [161117] . We note that repulsive short-range 3NF component 
is phenomenologically required to explain the properties of high density matter 3,7,8 
and its detailed information is one of the central question of interest. 

Since 3NF is originated by the fact that a nucleon is not a fundamental parti- 
cle, it is most desirable to determine 3NF directly from QCD, in particular at short 
distance where the dynamics of fundamental degrees of freedom (DoF), quarks and 
gluons, becomes essential. In this proceeding, we report our progress on first-principles 
calculations of 3NF using lattice QCD simulations 18,19 . Note that while there are 
lattice QCD works for three- (or more) baryon systems [20,21,22,23], they focus on 
the energies of the multi-baryon systems, and extracting 3NF is currently beyond their 
scope. 

As for the calculation of two-nucleon forces (2NF) from lattice QCD, an approach 
based on the NBS wave function has been proposed 24 , 25 , 26 , so that the potential is 
faithful to the phase shift by construction. Resultant (parity-even) 2NF are found to 
have desirable features such as attractive wells at long and medium distances and cen- 
tral repulsive cores at short distance. The method has been applied to general hadronic 
interactions 27,28,29,30,31,32,33,34,35 . The method itself has been also generalized 
to "time-dependent" HAL QCD method, so that the energy-independent (non-local) 
potential can be extracted without ground state saturation [36]. See Refs. [26,37, 
138] for recent reviews. In this report, we extend the method to three-nucleon (3N) 
systems, and perform the lattice QCD simulations for 3NF in the triton channel, 
[I,J ) — (1/2, l/2 + ) 18,19]. In particular, we update the results in Ref. [19] by 
employing the time-dependent HAL QCD method [5S] to suppress the systematic er- 
ror associated with excited states. We also utilize a novel contraction algorithm (unified 
contraction algorithm), which reduces the enormous computational cost for 3N corre- 
lators by a factor of 192 [39] . 

2 Formalism 



We first briefly explain the framework for the calculation of 2N potentials 26,37,38 . 
The central quantity in the HAL QCD method is the equal-time Nambu-Bethe-Salpeter 
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(NBS) wave function, 4>2N( r ) = (0\N(r)N(0)\E2N) i where \E 2 n) denotes the state of 
the 2N system at the energy of E 2 n m the center-of-mass frame, N the nucleon operator 
where spinor/fiavor indices are implicit. For simplicity, we consider the elastic region 
hereafter, while the method can be extended above inelastic threshold [31II4U] , The 
important property of the wave function is that it has a desirable asymptotic behavior, 

, sinffer - Ztt/2 + 5f ) , , 

02Aroc i — i ^, r=|r|^c», (1) 

where E 2 n — k 2 /(2fj.) with the reduced mass fi = mpf/2, I being the orbital angular 
momentum. Exploiting this feature, we define the (non-local) 2N potential, U 2 ff(r,r ), 
through the following Schrodinger equation, 

- 2^-<fev(r)+ / dr f/ 2 jv(r,r )(/> 2 jv(r ) = E 2N (f> 2N (r). (2) 

It is evident that, although U 2 n itself is not an observable, U 2 n is always faithful to 
the phase shift by construction. Another important property is that, while U 2 n could 
be energy-dependent in general, it was proven that one can construct U 2 n so that it 
becomes energy-independent 25,26,38,40 . 

In lattice QCD, the NBS wave function of the ground state, </>2JV> can ^ e extracted 
from the four-point correlator as 



G 2N (r,t-t Q ) = ^^(0|(iV(R + r)iV(R))(i) (N'N')(t )\0), (3) 

R 

— ^2A^2V~ B2Vt - to) , A° 2N = <E 2 %pW)|0>, (4) 

<f> 2N (r) = {0\N(r)N(0)\E%), (5) 

where E 2N denotes the energy of the ground state, N (N 1 ) the nucleon operator in the 
sink (source). In the practical lattice calculation, however, it is notoriously difficult to 
achieve the ground state saturation [3B]. This is because (i) the energy splitting between 
the ground state and excited states are getting smaller for larger lattice volume and (ii) 
signal to noise ratio (S/N) in lattice Monte Carlo simulation is ruined exponentially 
with increasing t. S/N also becomes exponentially worse for larger mass number in the 
system, and/or for lighter pion mass on the lattice. 

In order to overcome this problem, it is recently proposed to consider the 
time-dependent Schrodinger equation, 

- -^-ip2N{^,t) + / dr'U 2N (r,r')i> 2N (r,t) = -^V2jv(r,t), (6) 

where V^iV is imaginary-time NBS wave function defined by 

^ 2N (r,t) = G 2N (r,t)/e' 2mNt . (7) 

What is noteworthy is that, thanks to the energy-independence of U 2 n, Eq. © holds 
not only for the ground state but also for excited states simultaneously. Therefore, the 
ground state saturation is not required, which is a significant advantage of the potential 
approach to multi-baryon systems in lattice QCD. 

In the practical calculation, we perform the derivative expansion for the non-locality 
of the potential, U 2N (r,r') = [V c (r) + V T (r)S 12 + V LS {r)L ■ S + C(V 2 )] 5(r - r'), 
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where Vq, Vj- and Vls are the central, tensor and spin-orbit potentials, respectively. 
In Ref. [33] , the validity of this expansion is examined, and it is shown that the leading 
terms, Vc and Vt, dominate the potential at low energies. 

We now extend the method to the 3N system. We consider the imaginary-time NBS 
wave function of the 3N, ip3Ni T > Pi *)i defined by the six-point correlator as 

G 3N (r, p, t - t ) = i £)<0|(JV(xi)JV(x 2 )JV(x3))(t) (N'N'N>)(t )\0), (8) 



R 



^ 3N (r,p,t-t ) = G m (r,p,t-t )/e- 3m ^ t - t ^ (9) 

where R = (xi + x 2 + X3)/3, r = xi — x 2 , p = X3 — (xi + x 2 )/2 are the Jacobi 
coordinates. With the derivative expansion of the potentials, the NBS wave function 
can be converted to the potentials through the following Schrodinger equation, 



~y 2 r ~ + V V 2N (r ij ) + V 3NF (r,p) 



d 

^3N(r,P,t) = --^i>3N(r,P,t), (10) 



where V^jv( r »i) with r,j = Xj — Xj denotes the 2NF between (i, j')-pair, V37VF( r > p) the 
3NF, /i r = mjy/2, ^ip = 2m]\r/3 the reduced masses. If we calculate ip 3 ^(r, p,t), and 
if all V 2 7v( r y ) are obtained by (separate) lattice calculations for genuine 2N systems, 
we can extract Vsjyp (r,p) through Eq. (|10p . 

One of the greatest challenges in the study of multi-baryon systems on the lattice 
is that the computational cost of the correlators is exceptionally enormous. Actually, 
for larger mass number in the system, the cost for color/spinor contractions grows 
exponentially, and the cost for Wick contractions grows factorially: The total cost is the 
multiplication of both costs. In order to meet this challenge, we recently develop a novel 
algorithm in which color/spinor contractions and Wick contractions are considered 
simultaneously in a unified contraction index list [35]. In this "unified contraction 
algorithm," redundancies in the original contractions are eliminated, and a significant 
reduction in the computational cost is achieved, e.g., by a factor of 192 for 3 H and 
3 He correlators, and a factor of 20736 for the 4 He correlator. We note that, since a 
potential is independent of the choice of a nucleon operator at the source, we here 
employ the non-relativistic limit operator at the source to maximize the gain by the 
unified contraction algorithm. For the nucleon operator at the sink, which defines the 
NBS wave function (and correspondingly, the "scheme" of the potential 26,37,38 ), 
we employ the standard nucleon operator, N stc i = c a bc{<la C-ysq^qc, for both of 2NF 
and 3NF, so that they are determined on the same footing. 

In our first exploratory study of 3NF, we restrict the geometry of the 3N. More 
specifically, we consider the "linear setup" with p = 0, with which 3N are aligned 
linearly with equal spacings of r 2 = \r\/2. In this setup, the third nucleon is attached 
to (1, 2)-nucleon pair with only S-wave. Considering the total 3N quantum numbers of 
(/, J p ) = (1/2, l/2 + ), the triton channel, the wave function can be completely spanned 
by only three bases, which can be labeled by the quantum numbers of (l,2)-pair as 
1 Sq, 3 5i, 3 -Di. Therefore, the Schrodinger equation leads to the 3x3 coupled channel 
equations with the bases of ipi g , tp3g , ip3£ >1 . The reduction of the dimension of bases 
is expected to improve the S/N as well. It is worth mentioning that considering the 
linear setup is not an approximation: Among various geometric components of the 
wave function in the triton channel, we calculate the (exact) linear setup component 
as a convenient choice to study 3NF. While we can access only a part of 3NF from it, 
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we plan to extend the calculation to more general geometries step by step, toward the 
complete determination of the full 3NF. 

We consider the identification of genuine 3NF. It is a nontrivial work: Although 
both of parity-even and parity-odd 2NF are required to subtract 2NF part in Eq. (|10|) . 
parity-odd 2NF have not been obtained yet in lattice QCD. (See, however, our recent 
progress [35].) In order to resolve this issue, we consider the following channel, 

i>s = [ - Pt n t n i + Pt n i n t ~ ri t r HPt + n i n tPt + n tPt n i - n iPt n t] > ( n ) 

which is anti-symmetric in spin/isospin spaces for any 2N-pair. Combined with the 
Pauli-principle, it is automatically guaranteed that any 2N-pair couples with even 
parity only. Therefore, we can extract 3NF unambiguously using only parity-even 2NF. 
Note that no assumption on the choice of 3D-configuration of r, p is imposed in this 
argument, and we thus can take advantage of this feature for future 3NF calculations 
with various setup of 3D-geometries. 

3 Lattice QCD setup and Numerical results 

We employ Nf = 2 dynamical configurations with mean field improved clover fermion 
and RG-improved gauge action generated by CP-PACS Collaboration |41] , We use 
598 configurations at (3 = 1.95 and the lattice spacing of a~ 1 = 1.269(14) GeV, and 
the lattice size of V = I? x T = 16 3 x 32 corresponds to (2.5 fm) 3 box in physical 
spacial size. For u, d quark masses, we take the hopping parameter at the unitary 
point as k uc i — 0.13750, which corresponds to = 1.13 GeV, m^r = 2.15 GeV and 
= 2.31 GeV. Lattice simulations are performed at eleven physical points of the 
distance r-2- We use the wall quark source with Coulomb gauge fixing. In order to 
enhance the statistics, we perform the measurement on 32 wall sources using different 
time slices, and the forward and backward propagations are averaged. The results from 
both of total angular momentum J z = ±1/2 are averaged as well. For the sink time, 
we calculate for a wide range of 5 < (t — to) /a < 10. We evaluate Eq. (|10[1 at t being 
integer or half-integer, and we adopt the symmetric difference on the lattice for the time 
derivative. In the case of t of integer, ^^]\f(r,p,t) is obtained directly on the lattice, 
while ^■i/ , 3At(t', p, i) is obtained from ip^^(r,p,ti: 1). In the case of t of half-integer, 
both of ip3N( r , Pi t) an d ■§{' l l ) 3N( r j P^) are evaluated from ^>3jv(f, p, i ± 1/2). 

In Fig.[TJ we plot the radial part of each wave function of ipg — (—tpi g +ip3 g 1 ) /\/2, 
ipM = (^is + ifj&g^ / \/2 and tp3jj 1 obtained at (t — to) /a = 8. Here, we normalize 
the wave functions by the central value of ipsi r 2 = 0). What is noteworthy is that the 
wave functions are obtained with good precision, which is quite nontrivial for the 3N 
system. We observe that ipg overwhelms other wave functions. This indicates higher 
partial wave components are strongly suppressed, and the effect of the next leading 
order in the derivative expansion, spin-orbit forces, is suppressed in this lattice setup. 

We determine 3NF by subtracting 2NF from total potentials in the 3N system. 
Since we have only one channel (Eq. (|ll|l ) free from parity-odd 2NF, we can determine 
one type of 3NF. In this work, 3NF are effectively represented in a scalar-isoscalar 
functional form, which is often employed for phenomenological short-range 3NF [13] . 

In Fig. [2] we plot the preliminary results for the effective scalar-isoscalar 3NF at 
(t — to) /a — 7.5, 8.0, 8.5. The results from different sink times are consistent with each 
other within statistical fluctuations. This indicates that contaminations from excited 
states above inelastic threshold and/or the higher order terms in the velocity expansion 
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Fig. 1 3N wave functions in the triton chan- 
nel at (t — to) /a = 8. Circle (red), triangle 
(blue), square (green) points denote tjjg, ipMi 
i/>3£> , respectively. T2 is the distance between 
the center and edge in the linear setup. 



<t-t )/a=7.5 
(t-t )/a=8.0 i 
(t-t )/a=8.5 i 



0.5 1 

r 2 [fm] 

Fig. 2 The effective scalar-isoscalar 3NF 
in the triton channel with the linear setup. 
Green, red and blue points (with offset for 
visibility) are obtained at (t — to) /a = 7.5, 
8.0 and 8.5, respectively. 



of the potential are marginal. Fig. [2] shows that 3NF are small at the long distance 
region of r%. This is in accordance with the suppression of 27rE-3NF by the heavy pion. 
At the short distance region, however, an indication of repulsive 3NF is observed. We 
recall that a repulsive short-range 3NF is phenomenologically required to explain the 
properties of high density matter. Since multi-meson exchanges are strongly suppressed 
by the large quark mass, the origin of this short-range 3NF may be attributed to the 
quark and gluon dynamics directly. In fact, we recall that the short-range repulsive (or 
attractive) cores in the generalized two-baryon potentials are systematically calculated 
in lattice QCD in the flavor SU(3) limit, and the results are found to be well explained 
from the viewpoint of the Pauli exclusion principle in the quark level 28,34 . In this 
context, it is intuitive to expect that the 3N system is subject to extra Pauli repulsion 
effect, which could be an origin of the observed short-range repulsive 3NF. Further 
investigation along this line is certainly an interesting subject in future. It is also of 
interest that the analyses with operator product expansion [42] show that 3NF has a 
repulsive core at short distance. 

Evaluation of systematic errors are in progress. In particular, the effect of the 
discretization error is important to be investigated, since the nontrivial results are 
obtained at short distance. For this purpose, an explicit lattice simulation with a finer 
lattice is currently underway. Quark mass dependence of 3NF is certainly an important 
issue as well, since the lattice simulations are carried out only at single large quark 
mass. In the case of 2NF, short-range cores have the enhanced strength and broaden 
range by decreasing the quark mass [26]. We, therefore, would expect a significant quark 
mass dependence exist in short-range 3NF as well. In addition, long-range 27rE-3NF 
will emerge at lighter quark masses, in particular, at the physical point. Quantitative 
investigation by lattice simulations with lighter quark masses are currently underway. 



4 Conclusions and Outlook 



We have explored three- nucleon forces (3NF) in lattice QCD, utilizing the imaginary- 
time Nambu-Bethe-Salpeter wave function. The enormous computational cost is dras- 
tically reduced by the newly- developed "unified contraction algorithm." Using Nt = 2 
dynamical clover fermion configurations at a — 0.156 fm, V = 16 x 32 and = 1.13 
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GeV, we have studied 3NF in which three nucleons are aligned linearly with an equal 
spacing. Repulsive 3NF have been found at short distance in the triton channel. 

Currently, nuclear forces in lattice QCD are determined at rather heavy quark 
masses, which is considered to be a largest source of systematic errors. However, thanks 
to the significant theoretical and hardware development, it becomes possible to perform 
lattice simulations at the physical quark mass point. Generation of physical point gauge 
configurations with large volume has started, and nuclear forces will be subsequently 
studied [35]. While there remain various challenges, it is becoming within reach to de- 
termine realistic nuclear forces including few-baryon forces from first-principles lattice 
simulations, which will play an ultimate role in nuclear physics and astrophysics. 
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